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Abstract 


We offer a method for solving the fractional optimal control problems of 
multi-dimensional. We obtain a fractional derivative and multiplication op- 
erational matrix for Mott polynomials (M-polynomials). In the proposed 
method, the Caputo sense of the fractional derivative is applied on dynam- 
ical system. The main feature of this method is to reduce the problem into 
a system of algebraic equations in order to simplify it. We also show that 
by increasing the approximation points, the responses converge to the real 
answer. When the degree of fractional derivative approaches to 1, then the 
obtained solution approaches to the classical solution as well. 
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1 Introduction 


Over the last years, many researchers have developed fractional calculus and 
its applications in physics, chemistry, engineering, and so on; see [2, 7, 8, 
12]. A fractional dynamic system (FDs) is a system whose dynamics are 
expressed by fractional differential equations, and a fractional optimal control 
problem (FOCP) is an optimal control problem for an FDs [4, 26]. The 
theory of optimal control is an area in mathematics that has been developing 
for years, but the theory of fractional optimal control is a new area. Several 
authors studied this field [1, 3, 9, 22]. Some applications of FOCPs are given 
n [10, 19]. In most papers, one-dimensional FOCPs are considered where 
the problem is only with one state, one control variable, and one fractional 
differential equation [14, 15, 28]. In recent years, the use of different methods 
for solving optimal control problems has been considered by some researchers. 
A new numerical approach for solving FOCPs, including state and control 
inequality constraints using new biorthogonal multiwavelets, was made by 
Ashpazzadeh, Lakestani, and Yildirim [6]. In this paper, we consider FOCPs 
in the Caputo sense with multi-dimensional variables. In this method, the 
state and the control vectors do not necessarily have only one dimension. 
An FOCP can be defined with respect to different definitions of fractional 
derivatives. Most important types of fractional derivatives are the Riemann— 
Liouville and the Caputo. In this paper, we consider the multi-dimensional 
FOCP as follows: 
Determinate states X(t) € R” and controls U(t) € R*, which 


Min J(t, X(t n= f f(t, X(t), U(#))dt, (1) 
subject to 
D°a;(t) = gi(t, X(t), U(t)), 0<a<10<t<1li=1,...,n, (2) 
and satisfy the initial condition 
X(0) = Xo, (3) 
where 
X(t) = [x1(t),..-,n(t)?, Ut) = [un (t),--- ue)’, (4) 
Xo = [0,1 (t),---,2o,n(t)]7, and also f, gi : [0,1] x R” x R* > R are poly- 


nomial functions. The above problem reduces to a standard optimal control 
problem, when a = 1. 


The main purpose is to generalize the Mott operational matrix to frac- 
tional calculus. In the present work, we use Mott polynomials (MPs) for 
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solving FOCPs. The method consists of expanding the solution by MPs with 
the unknown coefficients. The properties of MPs are used to evaluate the 
unknown coefficients and find an approximate solution for X(t) and U(t) 
in the problem. Also, illustrative examples are included to demonstrate the 
applicability of the new approach. 

This article is organized as follows. Section 2 presents some of the pre- 
liminaries in fractional calculus. Section 3 describes the MPs and function 
approximation. We make an operational matrix for fractional integration, 
derivative, and multiplication by MPs in Section 4. In Section 5, we apply 
the MPs to solve multi-dimensional FOCPs. In Section 6, the convergence 
of the proposed method is discussed. In Section 7, the numerical examples 
are simulated to demonstrate the high performance of the proposed method. 
Finally, Section 8 concludes our work. 


2 Some preliminaries in fractional calculus 


This section provides some basic definitions of the fractional calculus. 


Definition 1. We define from [11] C,, = {f(t) : ne ] a t>0,f() = 
Pe — ) € C,}, where 
née Nand weR. 


Definition 2. The Riemann—Liouville fractional integral operator of order 
a> 0 of a function f € Cy, w > 1, is defined as [11] 


— s)*- il! oy 
eft) <4 Flay be) felds, a> 0.t>0, 
f(t), a=0, 


(5) 


and forn-L<a<n,neé€N,t>0,f € C™,, the fractional derivative of 
f(t) in the Caputo sense is defined as 


oDi f(t) = 1" D" f(t) 
: f(t - ayo F(s)ds, n-l<a<n,neN, 
{MW), a=nneN, 
(6) 
where D” is the nth-order derivative. 


Definition 3. For f € C,,u>—1,a, 8 > 0, we have [27] 


eG ale 1) pyte (7) 


[et) = ; 
a T(yta+1) 
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and forn-l<a<nneéN,t>0,f € Ci,u = —1, also, we have the 
following properties: 


1. pD°r* f(t) = f(t), 


ak 


n-1 
2. 1% 6D°f() = FO) -— YF OF) T, 
j=0 k} 


3,308 fia IeP Desa). 


3 MPs 


The MP s,,(t) is defined by (see [21, 23, 25]) 


#\" h($) 472k 
si) =(-1(5) o-oo ea 


where 


if n is even, 


1 
ae if n is odd. 
From the above definition, it is obvious that so(t) = 14 0. For this reason 
S,(t) is a Sheffer set [29]. The first few MPs are 


so(t) =1, 
1 
s1(t) = 2 
82(t) = qo 
s3(t) = oi a 
sa(t) = =e + att 
s(t) = - ~ ee : eS 


Now we define the vector A,, for r = 0,1,...,n. If r is odd, then for i = 
O, Leys 
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A, = bi,, (9) 
_ f0, i=21,le NU {0}, 
Bin = i otherwise, 70) 
=i” .)(=1) 
Ci, = (r—1)! (JCD ' 
' 2 (r — 2k; — 1)! 
ky = B;, By —1,.+.,0, 
r—1 r—-—1 
Bag a. PSO ews (11) 
If r is even, then for 7 = 0,1,...,n, 
A, = 6, , (12) 
, _fO i=21+1le NU {0}, (13) 
‘r | ci, otherwise, 
= (ZY ¢- 1 WO 
e 2 “(r— 2k —1)P 
ki, = 85,85 —1,-.-,0, 
st we ae a 
Bp = 5 qd 9 0,1,...,5. (14) 


Hence, $,(t) = A,T,(t) for (r = 0,1,...,n), here T(t) = [1,¢,...,é7]*. 
Then, we define an (r + 1) x (r + 1) lower triangular matrix A such that 
A=[Ao, A1,...,Ar]? and A;(i = 0,1,...,n) is a row vector of order r. 
As a result, 
Qn(t) = AT, (0), (15) 
where 
T 
Yn(t) = [so(t), s1(t),..., Sn(t)] : (16) 


3.1 Function approximation 


We recall here a theorem that was stated and proved in [20]. Suppose that 
H = L?(0, 1] is a Hilbert space and that {s0, s1,...,8n} is the MPs of degree 
n on the interval [0,1]. We define Y = Span{so,s51,...,5n}. Let f be an 
arbitrary element in H. Since Y is a finite-dimensional subspace of the space 
H, the function f has the best unique approximation on Y like f, € Y, that 
is, there exists f, € Y such that 
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If — fall2 < If — ylle, for all y € Y, (17) 
where ||f|l2 = \/(f, f) and (-,-) denotes the inner product. Since f, € Y, 
therefore f,, is a linear combination of the spanning basis of Y; that is, there 
exist n + 1 coefficients 


C =[co,¢1,---,Cn] ER (18) 


such that 


{O20 =) e50=C' o.@: (19) 
j=0 


where 


lf — fnlle > min. 


(20) 
Then C can be obtained by 
C= QO F@), Yni(t)), (21) 
where 
Q= (onl oalt)) =f ealOen Oat (22) 


Theorem 1. [20] Let X be an inner product space and let M #4 @ be a 
convex subset that is complete in the metric induced by the inner product. 
Then for every given x € X, there exists unique y € M such that 


§ = inf lea = lle — all (23) 


4 Mott operational matrices 


4.1 Mott operational matrix of the fractional integration 


In this section, we describe the MPs operational matrix of fractional integra- 
tion of the vector y,. The operational matrix can be approximated as 


Ign (t) ~ Pen(t), (24) 
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where P is the (n+ 1) x (n+ 1) Riemann-Liouville fractional operational 
matrix of integration. We construct P as follows: 


wazhue h(3) j a es 
I%s,(t) = (—1) (5) (i—1)! » (j) 1) =a 
i h(s i . 
2 rae (i —2k— 1)! Ti-2k+a+4+1) 
(25) 
Now we approximate t’~?*+° by n+ 1 terms of the Mott. basis 
EP RTs ST Balt), (26) 
j=0 
where 
b= Oe ye AE arta) de 
0 
h($) rg 
oe eae (7) (-))” 1 
xn (5) GU! Oe GID) (eae SOL 
(27) 
Therefore we have 
I%s,(t) ~ > Bis;(8), (28) 
j=0 
where 
—f4\%+9 
By =(-1)9 (5) @- 1G -1) 
h(S)A(S) vy. i 
(Cem 
Ee EEND) AR G—2k Ij SL 1 
1 
‘Ppa = Oh ae (72) 
and 
1 
Q; = (siteas(0) = fF assy Pat (30) 
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Finally, we obtain 


Boi Bee 


where P is called the MPs operational matrix of fractional integration. 


4.2 Mott operational matrix of the fractional derivative 


In this section, we describe the MPs operational matrix of fractional deriva- 
tive of the vector y,(t). The operational matrix can be approximated as 


oD On(t) = Ten(t), (32) 


where T is the (n+1) x (n+1) operational matrix of the fractional derivative. 
We construct T as follows: 


6 Dey, (t) = Woe | (¢—s)r03 = 9 (s)ds 


{POH 2 pl), 


T'(n— a) 
where * denotes the convolution product and 
eh (t) = D™ yn(t). 


Therefore we have 


a 1 =O — nm 
D“gn(t) = Tiga) * * D' yy (t) 

pe yr 

= ot y(t) = ——_t"- |! & AT (t 

T(n— a) * en(t) T'(n — a) : (t) 

Ae 1 — = n 

= ae 2) ee 
AD" 

[hw a C16 a oe 
n-@ 

= AD" |D"1, D*t,..., Dt") 

= Ap” 0! t& i fi-@ n} fro 

7 T(l-a) ’?TQ@-a) "" Tin-a+1) 

= AKD"T,, 
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where K is an (n+ 1) x (n+ 1) matrix where the entries are given by 


i! 
eo naeiea.. = 
0, t# J, 


i,j =0,1,...,n, 


and 
Tile Oe a 
Now we approximate t”~% by n+ 1 terms of the Mott basis 


Pe eee oye), 


where 
Bag fe *on( nar = Qn [em eary(tat 
oa [ POP ast lat 
0 
1 1 1 e 
sl 
Q | atl’n-at+2’ — ; 
and 


_f AT, (t)Tn(t)? AT dt = A [mW t)T,dtA’ = AHA’, 
where H is the well-known Hilbert matrix 

1 1 1 

, 3 3 i 
n 

1 1 1 { 

H= 2 3 4 n+2 

a Coe 


ntln+2n+1. Intl 
Finally, we obtain 
ODE Yn(t) ~ Ten(t), 
where 


T =AKD"E. (33) 
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Moreover, T is called the MPs operational matrix of fractional derivative. 


4.3 The operational matrix of multiplication 


The following property of the product of two Mott function vectors will be 
also used: 


Cw, (ie, (tb) go, C", (34) 


where C is the (n+1)x(n+1) multiplication operational matrix. To illustrate 
the calculation procedure, we let 


CT yn (t)en(t)? = [co(t), cx(t),... sealh)|. 


Now we approximate C7 yn(t)yn(t)” as follows: 


j=0 
where 
O; = G0, é,..-, Gin]. (36) 


Using (21), we obtain 


j=0 4=0 


where c¥ = (c;(t), ,(t)) and 


dix = (8;(t), x(t). (38) 
Therefore 
CT =C7D, (39) 
and 
C= |e ewe: (40) 


where D = [d,,] is a matrix of order (n +1) x (n+ 1) given by (38). Also, 
C? in (39) is given by 
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CF =cTp"}. (41) 
Hence, we get the operational matrix of multiplication as 


Co = [ees las pea (42) 


5 MPs for solving multi-dimensional FOCPs 


Using Theorem 1, we can approximate the state functions x;(t) and the 
control functions u,(t) as 


ni(t)© CT yalt), t= 1... (43) 


uj(t)* Brent), f= Leak, (44) 
where C;, B; € R°”*)*1. By using (32) and (43), we can write 


De Dep), telecom 


As a result, problem (1)—(4) reduces to 


1 
Minimize : F(t, CT Gn(t),---,Cn Pn(t), BT galt), ..-, BE Yn(t)) dt, 
0 
subject to 
CFT gn (t) = gi(t, CT On(t),.--,CZ Ont), Bi Gn(t),.--, BE Gn(t), 
with the initial conditions 
CF yn (0) = X05 _— dl svete WU (45) 


Since functions f,g; are polynomials, we get the following approximations: 


1 
| f(t, X(t), U(b))dt = F(Ci,...,Cn, Bi,..., Br), (46) 
0 
gilt, X(t), U(t)) oS G,(Ci, paral Cn, Br, missy Bx) en(t), t= 1, trey TN, 
(47) 
where F,G; : RO*+)xn x R(mt+1)xk _, R1x(m+1)_ For each i = 1,...,n, we 


can generate algebraic equations from (47) as follows: 
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1 
Gig = | (Orr = Gi(Ci, eens Bigc.ey Br) Pn(t)Byn(t))dt = 0, 
0 


and from (45), we set Gin = CP yn (0) — x04. 


Finally, the FOCP (1)-(4) has been reduced to a parameter optimization 
problem, which can be stated as follows: 


Find C; and B; that 


Min F(C,...,Cn, Bi,..., Br) 
S.t. Gog Cig Og Biases By) = 0, 
i=1,...,n,j =0,...,k. 


We define the Lagrange function for the above problem as 


L(C\,...,Cr, Bi, «+, Bry At, +++; An) 


k 
SF Chyie Cy Bigwvey Bn) YY ha gGeg (Ci axss Cas Big Be), 


—0 


ll 
an 
aS 


and consider the necessary conditions for the extremum and obtain the cor- 
responding system of algebraic equations 


OL OL 

= = — L=1 oe 
aC; 0, a; 0, a ’ MN, 
OL 
=, =0, j=l k 
OB; ? J ? ? 


These equations can be solved for C;,B;, and de by Newton’s iterative 
method. Then, we get the approximate value of the state functions «;(t) 
and the control functions u;(t) from (43) and (44), respectively. 


6 Error estimation and convergence analysis 


In the following theorem, the error estimation for the approximated functions 
will be expressed in terms of Gram determinant. For any given elements 
%1,02,..-.,%p in a Hilbert space H, the Gram determinant of these elements 
is defined as follows [20]: 
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ae ee Sos eee 
G(21,22,...,2n) = 1 si 2 - az n (48) 
(fn, £1) (Ln, LQ)... (Ln, Ln) 


Theorem 2. Suppose that H is a Hilbert space and that Y is a closed 
subspace of H such that dimY < oo and yi, y2,...,Yn is any basis for Y. Let 
x be an arbitrary element in H and let yo be the unique best approximation 
to x out of Y. Then from [20] 


G(x, y1, Ye, enh 5 Yn) 


Iz — yoll2 = ch Tee (49) 
where 
(z,@) (@,y1) --. (£,Yn) 
Ci Be (yt 2) (yay) a (vt Yn) (50) 


(i. Wien aed 


The following theorems illustrate that by increasing the number of MPs 
the error tends to zero. 


Theorem 3. Let f be an arbitrary element in H. Then the function f has 
the best unique approximation on Y like f,, € Y; that is, 


there exists f, € Y s.t. for ally¢Y: ||f — falle < |lf —ylle, (51) 
where ||f|l2 = /(f, f) and (-,-) denotes the inner product. Since fp, € Y, 
therefore f,, is a linear combination of the spanning basis of Y; that is, there 
are n+ 1 coefficients 


C = [co,C1,---, Cn] € R, (52) 


such that 


fO2 ho => 2 O=C' o.@, (53) 
j=0 


where 


lf — Frlle — min. 
(54) 


Then C can be obtained by 
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C = a ®), Yn(t)), 


where 


Q = (ynlt), Pat) = | ea 


Theorem 4. Suppose that f(t) € L?[0,1] is approximated by f,,(t) as 


fn(t) = > ciBi(t) = CT yn (t), 
i=0 
where C' and y,,(t) are defined respectively in (55) and (16). Then 
tim [If () — falO)llz2p0.y = 0. 
The error vector eI® of the operational matrix is given by 


el = lel? el? .0yel?)? = Pe) =1"e. ©. 


From (49) and Theorem 2, we have 


iakta_ Sop] — (GCE, sit), 81 (8), --- sult) 
' Sr] = ( G(so(t), 81(t),--- 8n(8)) ) 


Hence, according to (59) and (28), we have 


lleZll2 = |I%si(t) — D5 Bis, (t) 
j=0 


enc OlQ eo aeerrentces 


1 : n 
x pio 2k ta _ bsS5 
i+j-2L-2+a+1 28 
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(56) 


(57) 
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fp 0 =p et 


Caer) 
G(so(t), s1(£),---, 8n(t)) : 


(61) 


By considering Theorem 4 and using (61), one can conclude that by in- 
creasing the number of the Mott bases the vector ef® tends to zero. 


7 Numerical examples 


To demonstrate the applicability of the numerical scheme, we apply the 
present method for the following illustrative examples. 


Example 1. [16] Consider the following two-dimensional FOCP: 


Min Tu, 2] = 5 | (22(t) + w2(t))at, 


s.t. D°a(t) = —ax(t) + u(t), 


Our aim is to find the pair (#(t),u(t)), which minimizes the performance 
index J. The exact optimal solution of this problem for a = 1 is as follows: 
a(t) = cosh(V2t) + 6 sinh(V2t), 
u(t) = (1 + V28) cosh(/2t) + (V2 + 8) sinh(V28), 
cosh(V/2) + V2sinh(V2) 


a V2 cosh(V2) + sinh(V/2) | 


In Tables 1 and 2, the absolute errors of our proposed method and the 
method presented in [16] for a = 1 are given. We see that the absolute 
errors of our approach are less than that of the method presented in [16]. 
Moreover, in Table 3, we show the comparison of the approximate optimal 
value of objective functional for a = 1. It can be seen that the approximate 
value J tends to J* = 0.1929092980931791356. 
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Table 1: Comparison of the absolute errors of approximate optimal state for Example 
1 


t Present Method Method in [16] 


0.1 6.23e-9 2.11e-5 
0.2 4.60e-10 9.71e-6 
0.3 8.94e-9 4.08¢e-7 
0.4 1.01e-9 5.76¢-7 
0.5 9.33e-9 5.66e-6 
0.6 3.89e-10 9.25e-6 
0.7 8.67e-9 8.35e-6 
0.8 9.96e-10 4.34e-6 
0.9 5.83e-9 2.59-6 


Table 2: Comparison of the absolute errors of approximate optimal control for Example 
1 


t Present Method Method in [16] 


0.1 1.le-7 6.74e-6 
0.2 1.57e-7 3.17e-6 
0.3 3.49e-9 5.92e-7 
0.4 1.41e-7 7.10e-7 
0.5 3.84e-9 2.01e-6 
0.6 1.38¢e-7 2.71e-6 
0.7 9.80e-9 2.11e-6 
0.8 1.50e-7 8.59e-7 
0.9 1.08¢e-7 8.93-8 


Table 3: Comparison of the approximate optimal value of J for Example 1 


n Present Method Method in [16] 
J error of J J[16] error of J [16] 
7 0.1929092980931 7.6e-15 0.192909308499 1.04e-8 
8 0.1929092980931 7.6e-15 0.192909298458 3.64e-10 
9 0.1929092980931 4.5e-18 0.192909298107 1.41e-11 


Example 2. [13] Consider the following two-dimensional FOCP: 


Min Jlu,a] = a (a2 (t) + 3 + u?(t))dt, 


st. Da, (t) = —a1(t) + vo(t) + u(t), 
D“x2(t) _ —2x2(t), 
x1(0) = x9(0) =1. 


At a =1, the exact solutions are 
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3 —2t 
a1(t) = 0.018352eV + 2.48165e— V2" — a 
xo(t) = eae, 

2t 


u(t) = 0.044305" — 1.0279322e~ V2" + 7 


In Tables 4—6, the absolute errors when a = 1 and n = 3,4,6 are demon- 
strated, and in Figure 1 (a-c), the results for a = 0.9 and n = 3,4,6 are 
plotted. In Table 7, the comparison of our numerical results for the minimum 
values of J with different values of a at n = 6,7,8 with the results obtained 
in [24] is tabulated. In Figure 2 (a-c), we show that when a tends to 1, the 
approximate solutions tend to the exact solutions. In [13], a reasonable result 
was achieved with a large number of approximations (n = 64,128), while in 
the present work, we achieve a satisfactory result with at most six elements 
of the Mott basis, which demonstrates the efficiency of the new method even 
when the number of approximations is not so great. 


Table 4: Absolute errors of x1(t) for a = 1 at various choices of n for Example 2 


n 


t 3 4 6 

0.1 6.2546 x 10-4 7.4441 x 1074 1.5002 x 107° 
0.2 4.5464 x 1073 9.0972 x 10-4 8.3446 x 10-6 
0.3 3.9041 x 1073 1.4174 x 107° 1.0492 x 107° 
0.4 1.1462 x 10-3 7.5444 x 1074 5.1504 x 1076 
0.5 1.9425 x 10-3 8.7140 x 1074 1.2149 x 107-5 
0.6 4.0841 x 1073 3.1893 x 1074 1.0835 x 107° 
0.7 4.3836 x 1073 5.3794 x 10-4 1.2403 x 10-5 
0.8 2.2303 x 1073 1.0500 x 107 5.0792 x 10-6 
0.9 2.7809 x 10-8 3.6221 x 10-4 1.6075 x 10-5 


Table 5: Absolute errors of x2(t) for a = 1 at various 


n 


choices of n 


for Example 2 


t 3 4 6 

0.1 1.8817 x 10-4 1.0437 x 1073 1.3953 x 1075 
0.2 1.0761 x 10-2 1.4093 x 10-3 7.5730 x 10-6 
0.3 1.0648 x 107? 1.3573 x 1074 9.8425 x 10-6 
0.4 4.7251 x 1073 1.0583 x 1073 4.5687 x 10-6 
0.5 3.0146 x 1073 1.3381 x 1073 1.1291 x 107° 
0.6 9.3023 x 1073 5.9271 x 10-4 1.2293 x 10~® 
0.7 1.1462 x 107? 6.9627 x 10-4 1.1414 x 1075 
0.8 7.3019 x 10-3 1.5621 x 107-3 4.8958 x 10~® 
0.9 4.9714 x 1073 6.4096 x 1074 1.4835 x 107-5 
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Table 6: Absolute errors of u(t) for a = 1 at various choices of n for Example 2 


Alavi, Haghighi, Yari and Soltanian 


n 


t 3 4 6 
0.1 6.6070 x 10-4 2.3079 x 10-4 5.0375 x 10-6 
0.2 1.7968 x 10-4 3.1360 x 10-4 2.7248 x 10-6 
0.3 7.8105 x 10-4 3.1989 x 10-5 3.5467 x 10-6 
0.4 9.3073 x 10-4 2.3164 x 107+ 1.6400 x 107 
0.5 6.1527 x 10-4 2.9416 x 1074 4.0622 x 10-6 
0.6 1.9916 x 107° 1.3128 x 1074 4.4423 x 1077 
0.7 7.0047 x 10-4 1.5236 x 10-4 4.1084 x 10-6 
0.8 1.0445 x 1073 3.4459 x 1074 1.7637 x 10~® 
0.9 5.7713 x 10-4 1.4011 x 10-4 5.3528 x 10-6 
Table 7: Values of J when a approaches to 1, for Example 2 
PresentMethod Method of [29] 

a n=6 n=T7 n=8 n=6 n=7 n=8 
0.6 0.329 09 0.329 09 0.329 10 

0.7 0.351 62 0.351 62 0.35162 

0.8 0.376 27 0.376 27 0.376 27 0.378 33 0.377 56 0.377 17 
0.9 0.403 08 0.403 08 0.403 08 0.403 98 0.403 66 0.403 46 
0.99 0.429 00 0.429 00 0.429 00 0.429 09 0.429 05 0.429 04 
0.9999 0.43195 0.43195 0.43195 0.431 96 0.431 96 0.43196 
1 0.43198 0.43198 0.43198 0.43199 0.43198 0.43198 


Example 3. [1, 5, 22] Consider the following FOCP: 


Min Jlu,2] = 5 | (x?(t) + w2(t))at, 


st. D°ax(t) =ta(t)+ u(t), O<a<l, 
x(0) =1. 


This problem has the exact solution J = 0.4842676962287272. In Figure 
3 (a,b), the obtained results of the variables x(t) and u(t) are plotted for 
different values of a. In Figure 4 (a,b), the comparison between the exact 
solutions and the proposed method is plotted for n = 8 anda = 1. In 
Table 8, the comparison of our numerical results for the minimum values of 
J for a = 1,0.99 at n = 5 with the results obtained in [18, 24] is tabulated. 
Obviously, our estimated results are in good agreement with them. 
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Figure 1: The behavior of the approximate solutions of Example 2 for n = 3, 4,6 and 
a = 0.9, with exact solution. (a)x1(t), (b)x2(t), (c)u(t). 
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(c) 
Figure 2: The behavior of the approximate solutions of Example 2 for n = 6 and 
a = 0.9,0.99,1. (a)xi(t), (b)x2(t), (c)u(t). 
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Table 8: The estimated values of J for different values of a and n = 5 for Example 3 
a Present Method Method of [24] Method of [18] 


1 0.484267 0.484268 0.484268 
0.99 0.483461 0.483468 0.483463 
0.9 0.475874 0.476067 0.475883 


0 01 02 03 04 05 06 O7 O08 09 1 “o 01 02 03 04 05 06 07 08 09 1 
t t 


(a) (b) 


Figure 3: The behavior of the approximate solutions of Example 3 for n = 8 and 
a = 0.75, 0.85, 0.95, 1, with exact solution. (a)ax(t), (b)u(t). 
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Figure 4: Comparison between the exact solutions and the proposed method of Exam- 
ple 3 for n = 8 and a=1. (a)a(t), (b)u(t). 


Example 4. [30] In this example, the vibration of a spring-mass-damper sys- 
tem subjected to an external force is considered. In particular, the response 
to harmonic excitations, impulses, and step forcing functions is examined. In 
many environments, rotating machinery, motors, and so on cause periodic 
motions of structures to induce vibrations into other mechanical devices and 
structures nearby. On summing the forces, the equation for the forced vibra- 
tion of the system in Figure 5 is obtained. It is common to approximate the 
driving forces F(t) 
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mi(t) + ca#(t) + ka(t) = F(), 


where m,c, and k are fixed numbers. Also F(t) represents the control force 
derived from the action of an actuator force represented by F(t) = bu(t), 
where 0 is a fixed number. The linear regulator problem has a specific ap- 
plication in the vibration suppression, and the performance index for this 
problem is defined as 


Fit) 


te mg |—> Fit) 


Figure 5: (a) Schematic of the forced mass-damper system assuming no friction on 
the surface and (b) free body diagram of the system of part (a) for Example 4 


where to and ty are initial and final times, respectively. We introduce the 
usual state variable notation 


T= 2, 
L2 = § DF x(t). 


Then, the equation of motion in a order can be written as 


b 
U. 
m 
By selecting tp = 0, ts = 1, c= 2, andm =a=b=k =1, the boundary 
conditions for this example are considered as 


‘en = oDFx(t)(0), 
(1) = 9DPx(t)(1), 
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or 


For a = 1, we obtain the following exact solutions: 


1393 676 76 
L(t) =|j,-¢ ae Nein i67") 


e152" — ____e 152 
423 423 


16157 sr, 1134 _ sey 

UN) sag aa 
O57 er, 2591 _aer, 

152 LB? 


— ere" + 163° 


(cos E74), 


and Je = 13.00484741498823, where x;-(t) = x(t). In Table 9, the results 
of J for different values of a and n are listed. It is seen that with the increase 
in the number of the Mott basis, the approximate value of J converges to 
the exact solution. Also Tables 10-12 demonstrate the approximation of 
x1(t), v(t), and u(t) for different values of n with a = 1. Figure 6 (a-d) 
illustrates the behavior of state variables x(t), x2(t), control variable u(t) 
and performance index J, respectively, for n = 7 and for different values of 
a with the exact solutions. Table 13 demonstrates the approximation of J, 
for n = 8 and different values of a. Also in Table 14, the absolute errors of 
J, x(t), v(t), and u(t) for Example 4 are calculated at various choices of n 
and a= 1. 


Table 9: Approximate solutions of J for different values of a and n for Example 4 


a 
n 1 0.99 0.9 
3 13.01904761904762 12.644 81 10.531 56 
5 13.0048478586403 12.638 42 10.209 67 
7 13.00484741498875 12.649 19 10.210 82 
9 13.00484741498822 12.654 24 10.220 66 
11 13.00484741498823 12.656 99 10.227 68 
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Table 10: Absolute errors of 21(t) for a = 1 at various choices of n for Example 4 


n 


t 5 7 11 
0 0 0 0 

0.1 8.0624 x 10-8 2.8936 x 10-9 3.3240 x 10-8 
0.2 8.9073 x 10-6 3.8441 x 10-9 2.0139 x 10-8 
0.3 1.5140 x 10° 5.1600 x 10-9 5.0016 x 10-8 
0.4 1.1564 x 10-5 2.6349 x 107° 8.2645 x 10-18 
0.5 1.2533 x 1075 7.1014 x 107° 3.1086 x 10-15 
0.6 4.7572 x 10~§ 1.6704 x 10-9 8.1934 x 10-18 
0.7 4.5973 x 107° 5.2971 x 107° 4.8572 x 10718 
0.8 7.9685 x 1076 3.1233 x 107° 1.9929 x 10-18 
0.9 3.9680 x 10~% 2.7298 x 10-9 3.1830 x 1078 
1 1.3322 x 107-15 1.2212 x 10715 4.4409 x 10716 


Table 11: Absolute errors of x2(t) for a 


= 1 at various choices of n for Example 4 


n 


t 5 7 11 
0 0 0 0 

0.1 7.8816 x 10-5 2.6708 x 10-8 1.0362 x 10-1! 
0.2 6.2887 x 1075 6.7790 x 10-8 1.4781 x 1071! 
0.3 1.2314 x 1074 4.6618 x 1078 1.1712 x 10-4 
0.4 6.2434 x 107-5 8.3958 x 1078 5.3897 x 107! 
0.5 4.2027 x 10-5 7.2681 x 1079 1.6518 x 1071! 
0.6 1.003 x 10-4 8.4910 x 10-8 5.3266 x 10-¥? 
0.7 7.2908 x 10-5 3.2740 x 10-8 1.1480 x 10-44 
0.8 8.3314 x 10-8 6.6651 x 10-8 1.4331 x 107" 
0.9 5.8537 x 1075 1.7511 x 10-8 9.9403 x 10-1? 
1 7.7715 x 10716 1.2212 x 10715 7.7716 x 10716 


Table 12: Absolute errors of u(t) for a = 1 at various choices of n for Example 4 


n 


t 5 7 11 
0 3.9149 x 10-8 4.0395 x 10-6 1.2538 x 1079 
0.1 9.5758 x 10-3 5795 x 10-% 2.1342 x 10710 
0.2 1.3556 x 1073 6.0445 x 1077 2.8410 x 1071! 
0.3 1.7251 x 1074 .1567 x 1078 2.5282 x 10710 
0.4 8.5641 x 1074 2.3114 x 1077 2.8751 x 10710 
0.5 1.0077 x 10-3 .1500 x 107° 3.2260 x 10714 
0.6 3.5723 x 10-4 3.5077 x 10-7 3.0539 x 10719 
0.7 5.0337 x 10-4 .0067 x 10~® 2.0282 x 10719 
0.8 8.2804 x 10-4 6.34 x 10-7 8.3564 x 10714 
0.9 1.5161 x 10-4 3112 x 1076 2.4373 x 10710 
1 1.1831 x 1078 3.307 x 107 1.1912 x 107-9 
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Table 13: Approximation values of J for n = 8 and different values of a 


a=1 
n 
1 13.004 847 414 988 234 117 130 714 182 4 
0.99 12.676 065 068 414 313 965 482 218 5176 
0.9 10.416 927 546 051 989 730 195 819 221 4 
0.8 9.004 530 163 410 743 299 622 066 577 65 
0.7 8.363 433 733 210 117 712 725 797 350 86 
Table 14: Absolute errors of J,x21,22,u for a = 1 and different values of n 
array 
n error of J error of 2 error of x error of u 
11 7.6085 x 10716 4.846 x 10-18 1.04 x 10-4 2.792 x 19710 
9 7.5908 x 10-15 3.974 x 10-!? 7.744 x 10-" 1.9 x 10-° 
8 3.8474 x 10-15 2.729 x 10-10 4.563 x 10-9 9.6 x 10-8 
7 5.2381 x 10-12 3.908 x 10-° 5.603 x 10-8 1.024 x 10 
5 4.4365 x 10-7 7.474 x 10-6 7.234 x 10-5 9.419 x 10-4 


0 Of 02 03 04 05 06 O7 O08 09 Oo Of 02 03 04 05 06 O7 O08 o9 1 
t t 


01 02 038 04 O05 06 O07 O08 o9 1 
t 


(c) (d) 


Figure 6: The behavior of the approximate solutions of Example 4 for n = 7 and 
a = 0.88, 0.90, 0.99, 1, with exact solution. (a)x1(t), (b)x2(t), (c)u(t), (d)J. 
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8 Conclusions 


In this paper, a new numerical method has been derived to find the approxi- 
mate solutions of the multi-dimensional FOCPs; this numerical method uses 
MPs. The Mott fractional integration matrix reduced the FOCP into an 
equivalent integral problem. By using the Mott fractional derivative and 
multiplication matrix, we can transform the equivalent functional integral 
equation problem into an algebraic system of equations, where this new prob- 
lem is most easy to solve. Some examples are presented to demonstrate the 
validity and applicability of the new method. MATLAB (2018) is used for 
computations in this study. 
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